Exactly solvable phase oscillator models with synchronization dynamics 
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Populations of phase oscillators interacting globally through a general coupling function f(x) have 
been considered. In the absence of precessing frequencies and for odd-coupling functions there 
exists a Lyapunov functional and the probability density evolves toward stable stationary states 
described by an equilibrium measure. We have then proposed a family of exactly solvable models 
with singular couplings which synchronize more easily as the coupling becomes less singular. The 
stationary solutions of the least singular coupling considered, f(x) — sign(is), have been found 
analytically in terms of elliptic functions. This last case is one of the few non trivial models for 
synchronization dynamics which can be analytically solved. 



The dynamics of systems consisting of large sets of mu- 
tual interacting units is a very interesting problem in sta- 
tistical physics. In some cases, these units can be thought 
of as subsystems characterized by hidden internal degrees 
of freedom obeying their own dynamics. This is the case 
of large populations of individuals forming complex sys- 
tems of interest in interdisciplinary fields such as biol- 
ogy, economy, neurophysiology and ecology ^J^] . On the 
other hand, the subsystems can be seen as entities inter- 
acting with each other in the presence of quenched dis- 
ordered external fields. The effect of the external forces 
is to pump energy into the system thereby leading to 
nonlinear oscillatory behavior. This last case describes 
charged wave instabilities in plasmas dynamical re- 
sponse of Josephson junction arrays in the presence of an 
external ac field Q or nonlinear oscillations and coherent 
motion of magnetized domains in strongly coupled mag- 
netic systems ||. It was realized by Kuramoto Q that 
one could realize this rich dynamical behavior by con- 
sidering mean-field models of phase oscillators with ran- 
domly quenched natural frequencies. The crucial feature 
of the Kuramoto model is that the random precession 
frequencies play the role of external forces which pump 
energy into the system. This leads to dynamical insta- 
bilities and synchronization in low dissipation regimes. 
Moreover, the mean-field character of the model captures 
the essential mechanism which leads to synchronized dy- 
namics while retaining mathematical simplicity. 

In the simplest setting, we consider the dynamics of 
a system of nonlinearly globally coupled phase oscilla- 
tors with random frequencies u>i taken from a distribu- 
tion g{oj) and subject to external independent white noise 
sources r\i (of strength V2T) : 
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i = 1, . . . , N. Here <pi(i) denotes the ith oscillator phase, 
K > represents the coupling strength and / is a generic 
real function of periodicity 2tt. In Fourier space, the 
latter can be decomposed as f(x) = X)^L-oo a ne mx , with 

Important work on the Kuramoto model with a gener- 
alized coupling function was carried out by Daido |?J who 
introduced the concept of order function and Crawford 
H who proposed scaling behavior of bifurcating branches 
from the incoherent solution. These works introduced 
rather general theories and methods which, however, had 
inherent limitations (zero temperature for the order func- 
tion theory and vicinity to bifurcation points in Craw- 
ford's work). Thus it would be desirable to have exact 
results for solvable models where such theories could be 
checked and extended. 

The purpose of this letter is to show some exact re- 
sults on the synchronization dynamics of simple models 
of oscillators. Exact results are very difficult to obtain 
in models with general smooth coupling functions f(x). 
Here we will follow a different strategy by analyzing mod- 
els with singular couplings which can be mapped into well 
known physical problems. Later we will go further and 
consider models with smoother f(x) such as the Daido 
coupling which we analytically solve. While our analysis 
is restricted to models without disorder, the techniques 
developed here could be extended to disordered cases || . 

The moment approach and general equations. A sim- 
ple way to study the dynamics of the Kuramoto model 
with general coupling functions (0) is to use the recently 
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introduced moment approach (T^] and define, H™ = 
jfYlj=i < elfc< ^ J > ' wnere the brackets < ... > de- 



note average with respect to the external noise and the 
overbar denotes average with respect to the random os- 
cillator frequency. The equations of motion for the mo- 
ments read 



dt 
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The following differential equation for the generating 



function, g(x,y,t) = V fe= 
easy to derive from (0) 
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dx 2 dxdy 
Here v(x,t) is the drift velocity, defined by 

oo 

v(x,t)=-K J2 a » H -ne mx 
n=— oo 

= -X / /(z'M* - -t', 0, t) dx'. (4) 



It is a simple matter to check that this velocity satisfies 
v(x,t) = —K H(x — Q e t,t), where H(x,t) is Daido's 
order function |7) and fl e is the synchronization fre- 
quency (which is zero for stationary states). On the 
other hand, the generating function g(x, y, t) is related 
to the one-oscillator probability density p(x, u, t) through 
g{x, y, t) = e yuJ p(x, u, t) g(u) dw. We derive from (|) 
the usual nonlinear Fokker-Planck equation for p(x,u>,t) 
wherever g(u>) ^ 
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Existence of a Lyapunov functional. The solution of 
(||) is generally quite complicated: only in special cases 
it is possible to work out explicit results. Here we want 
to analyze under which conditions it is possible to find 
a Lyapunov functional for the dynamics. This is an im- 
portant result since such Lyapunov functional explicitly 
yields a functional equation for the stationary states as 
well as stability results. In particular we will find a Lya- 
punov functional when the coupling f(x) is an odd func- 
tion (i.e. if detailed balance is satisfied) and there is no 
frequency disorder [i.e. g(oj) — 5 (u>)]. Then there exists 
a partition function which characterizes the thermody- 
namic properties of the stationary states. 

Let us define the potential function V(x,t) = 
j_ v(s, t) ds and restrict ourselves to finding stationary 
solutions; rotating wave solutions may be reduced to this 
case after moving to a rotating frame. It is possible to 
show that the stationary solutions should have the form 



p(x, uj) = Z 1 e 



J 



j, I exp 



(x - s)lo + V(x) - V(s) 
T 



ds, (6) 



where Z and J are independent of x. We now im- 
pose the condition that p be a 27r-periodic function of 
x, use the symmetry properties of the drift and find 
the probability flux J as a function of Z: J/T — 
2Z- 1 sinh(27rT- 1 tj)/ e -^ +x ^ +v ^/ T dx. Z can be 
found from the normalization condition pdx = 1 . An 
interesting case corresponds to the case without disor- 
der, g(uj) = S(lo), for which J = 0, p(x, 0, t) = g(x, y, t) = 
g(x, t). We can obtain the stationary solutions (subscript 
zero below) by solving the system of equations: 



g (x) = Z 1 e Ys ^ L 



Z = 



e~^~ ds . 



v (x) 



f(x - x')g (x')dx' . 
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It is easy to check that, for this type of potential solu- 
tions, there is a Lyapunov functional [ |li"|] defined by the 
relative entropy, 



v(t) = 



g(x,t) In - 



g{x,t) 



g(x,t) 

V(x,t)- M (t 
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, . dV(x,t) , 
g(x,t) — — — dx. 
dt 
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(10) 
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Direct calculation (to be presented elsewhere) shows that 
rj{t) is bounded from below, rj'(t) < and that g tends 
to a stationary solution which is proportional to g. We 
then conclude that, for odd coupling functions and in ab- 
sence of disordered frequencies, the stationary state (Q) 
corresponds to the thermodynamic equilibrium state of 
a system of N oscillators interacting through the Hamil- 
tonian 



N ^ 



(12) 



where e(x) is a pair interaction energy 27r-periodic func- 
tion defined by e{x) = f f(s)ds. Note that V(x) — 
—K J_ g(x)e(x)dx + constant. Thus a computation 
of the partition function for ([TJ]) yields the stationary 
states. It is important to stress that potential solutions 
of the type (fjj) are no longer stationary solutions of the 
dynamics in the presence of disordered frequencies uii |ic| ] 
or when detailed balance is violated (i.e. f(x) is not an 
odd function) . The physical meaning of these two condi- 
tions is quite clear: suppose a single oscillator performcs 
a global rotation of angle 2tt. The global energy of the 
system, given by jl2]), does not change because e(x) is 
27r-periodic. Thus the coupling strength does not exert 
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work into the system of oscillators when a rotation path 
(<t>i — > 4>i + 2ir) is performed. Indeed, the amount of work 
is W = J^ l+2lT [LUi — J2f=i f( x ~^j)]^ x which vanishes 
only if LUi = and f(x) = —/(—a;). 

Singular coupling models. For general coupling func- 
tions, with infinitely many non-vanishing Fourier modes 
a n , it is difficult to find explicit analytical expressions 
for the stationary states. But calculations turn out to be 
much simpler in the case of singular coupling functions. 
Here we analyze three different cases of singular interac- 
tions: a) f(x) = S'(x); b) f(x) — S(x) and c) f(x) — 
sign(x). In the first case we have v(x) — —Kdg/dx, 



dg 

dt 



K 



d / dg 
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dx 2 
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This equation is related to porous media systems with 
the additional 27r-periodicity condition for g. It is easy 
to check that rj(t) — (1/2) J_ g 2 {x,t)dx is a Lyapunov 
functional. In fact, it is bounded from below and satis- 
fies, rf{t) = - J^(T + K Q g) (dg/dx) 2 dx < 0. Then the 
incoherent solution g(x) = l/2w is globally stable and 
model a) does not syncronize at any T. 

Case b) is described by the Burgers equation 



dg 
dt 



2K g 



dx 



T 



dx 2 



(14) 



for a 27r-periodic g, and it can be solved by the Hopf- 
Cole transformation. rj(t) defined for case a) is also 
a Lyapunov functional for case b), but now rf (t) = 
— T f_ 7r (dg/dx) 2 dx < 0. At finite T incoherence is again 
stable. But at zero temperature the oscillators remain 
synchronized forever if they are initially synchronized. 

Analytical solution of the model with Daido coupling. 
Model c) is analyzed carefully. This model was originally 
proposed by Daido in the general disordered case but we 
will solve it for g{u>) = S(lo). Since f(x) is odd, it is in 
principle possible to find the stationary states by solving 
the functional equations (§)-(§)■ Such a calculation is 
quite involved and here we follow a different and novel 
approach. The dynamical equations for model c) are non- 
local because the drift velocity (0) is 



v(x, t) = -K / [g{x -£,t)-g(x-Z + tt, (15) 
Jo 

but they become local after inserting the definitions 



p(x, t) = 2K [g(x,t)-g(x + ir, t)] , 
a{x,t) = 2K a [g[x, t) + g(x + tt, t)]. 

The new system is 

d(av) 



dp 
dt 

da d t dv\ 
-dt = ^ iv ] 
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The drift velocity v(x,t) is given by dv/dx = —p. All 
functions are 27r-periodic and it is easy to check that 
both p and v are antiperiodic if x — * x + tt, while a is 
7r-periodic. The stationary solutions of ([n])-([l8j) satisfy 
the following equations 



T 



dx 2 
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where L is a positive constant, 
once yielding 



L (19) 
0, (20) 
may be integrated 



1 / dw\ w 2 w 

2 [it) + T"T = ' 



(21) 



where w(£) = v(x)/y/L,x — T^/^/Z and £ is a new 
constant. ( f2l| ) can be easily solved by quadratures in 
terms of Jacobi elliptic functions p^ |. The periodicity 
properties of the drift v yields finitely many solutions 
(characterized by their value of £ and an odd integer n) 
for a fixed value of Ko/T. The relation between £ , n and 
K /T is: 



An 2 K(m) 



E{m) 



VI - 2£K{m) 
1 + s/l - 2£ 



= f- (22) 



Here K(m) and E(m) are the complete elliptic inte- 
grals of the first and second kind with parameter m — 
( 1 - \/l — 2£) / ( 1 + VI - 2£) , respectively. For every ad- 
missible value of £ and n, with Ko/T fixed, the drift 
velocity and the probability density can be found from 
@, ©, ©, © and the condition /„* o{x)dx = 2K Q . 
They are 



AnT^/mKim) 
v[x) = sou , 



K(m) 
2tt 



l + Vl-2£ 



to sn u 



to cnu aim 



E(m) 



Vl-2£ K(m) 



(23) 



(24) 



Here u = 2nK{m) (x — x )/tt, xq = constant, while sn, cn 
and dn are Jacobi elliptic functions defined in The 
details of these calculations will be presented elsewhere. 

As Ko/T — > oo the number of possible stationary 
solutions becomes infinite. Each different stationary 
state belongs to a different solution branch which bi- 
furcates from incoherence at a critical value of the cou- 
pling Kq = n 2 T7r/2, where n — 2p — 1 and p > 1 is an 
integer, see figure 1. The first branch bifurcates from 
incoherence at Ko/T = tt/2 and remains stable for all 
larger K /T. In the limit K /T -> oo, v(x)/K = 
sign(x — xo), g{x) = S(x — xq — tt), and full synchro- 
nization is achieved. A convenient synchronization order 
parameter r is defined through the global magnetization 
M 
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FIG. 1. Bifurcation diagram £ versus 2Ko/(ttT) for sta- 
tionary synchronized solutions branching off from incoherence 
£ — at the square of the odd integer numbers. 
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FIG. 2. Binder parameter g as a function of T for dif- 
ferent sizes (in the low T region, from bottom to top) 
N — 50, 100, 500, 1000. The curves become steeper for larger 
values of N. The theoretical value of the bifurcation temper- 
ature T — 2 /it is marked in the figure. 

The order parameter can be calculated from Monte 
Carlo simulations of the Hamiltonian ( [l2| ) with e(x) = \x\ 
and Kq ~ 1. The simulations use the heat-bath algo- 
rithm with a random sequential updating of the phases. 
The transition temperature corresponding to the first 
branch is easily obtained through standard finite-size 
scaling methods. Consider the kurtosis (or Binder pa- 
rameter) for the synchronization parameter r, g = i(3 — 



)i where < ... > is the standard configurational 
average [weighted with the usual Boltzmann-Gibbs fac- 
tor, exp(—f3H), and H is given by (U|)]. The curves for g 
are shown in figure 2 for different sizes. Note that these 
curves (specially for N = 50, 100, 500, data for N = 1000 
is more noisy) intersect at a common point characterizing 
the bifurcation temperature. 

Summarizing, we have considered models of oscilla- 
tors interacting through a general coupling function f(x). 
In the absence of precessing frequencies and for odd- 
coupling functions there exists a Lyapunov functional. 
Then the probability density evolves toward stable sta- 
tionary states which can be described by an equilibrium 
measure. In these cases the coupling force does not ex- 
ert work into the system and the dynamics is purely re- 
laxational. We have then proposed a family of exactly 
solvable models with singular couplings. The oscilla- 
tors may synchronize more easily for less singular cou- 
plings. In particular, case a), f(x) = S'(x), never syn- 
chronizes, case b), f(x) = S(x), retains synchronization 
only at zero temperature and case c), f(x) = sign(x), 
syncronizes at a finite temperature. This last case is 
one of the few non trivial models for synchronization 
dynamics which can be analytically solved. Stationary 
synchronized phases bifurcate from incoherence at a fi- 
nite temperature which corresponds to a thermodynamic 
second order phase transition. It would be very interest- 
ing to extend all our considerations in the presence of u>i 
for models with synchronization transition at finite tem- 
perature where we expect the appearance of dynamical 
instabilities as well as nonlinear behavior. 
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